Many-body coherent destruction of tunneling in photonic lattices 
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An optical realization of the phenomenon of many-body coherent destruction of tunneling, recently 
predicted for interacting many-boson systems by Gong, Molina and Hanggi [Phys. Rev. Lett. 103, 
133002 (2009)], is proposed for light transport in engineered waveguide arrays. The optical system 
enables a direct visualization in Fock space of the many-body tunneling control process. 
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The phenomena of coherent destruction of tunneling 
(CDT) [1] and dynamic localization (DL) [2] represent 
seminal results in studies of quantum dynamical con- 
trol. CDT and DL have been actively studied in sev- 
eral fields of physics [3], culminating in recent experi- 
ments on strongly driven Bose-Einstein condensates in 
optical lattices [U |S] which provided a very direct obser- 
vation of single-particle CDT and DL for matter waves. 
Many-body generalizations of CDT have attracted con- 
siderable interest in recent years, and experimentally |6l- 
111] . CDT for many interacting particles in a double well 
and its interplay with self trapping has been discussed 
in [5], whereas shaking- induced renormalization of the 
tunneling parameter and transition between superfluid 
and insulator states was predicted and observed in Refs. 
[71III]- In Ref.0, Gong, Molina and Hanggi recently pro- 
posed a novel route to CDT by considering a monochro- 
matic fast modulation of the self- interaction strength of a 
Bose-Einstein condensate in a double well potential in the 
framework of a two-site Bose-Hubbard model. Very inter- 
estingly, the modulation can be tuned such that only an 
arbitrarily, a priori prescribed number of particles are al- 
lowed to tunnel. In another physical context, the idea of 
CDT was proposed [TH [O] and experimentally observed 
[l4l [15] for light waves in coupled-waveguide structures, 
in which spatial propagation of light mimics the tempo- 
ral dynamics of a a driven quantum particle in a bistable 
potential [16j . The optical analogues studied in previous 
works [T2HT5] refer mostly to single-particle CDT. In the 
optical directional coupler proposed in Ref.[T3], the ad- 
dition of a Kerr nonlinearity enables to mimic CDT in 
presence of particle interaction in the mean-field limit, 
however such a scheme is not capable of simulating the 
two-mode Bose-Hubbard model. In a recent work [17] , it 
was shown that a bosonic junction, originally realized us- 
ing Bose-Einstein condensates in a double well potential 
[18j . can be simulated using engineered waveguide lat- 
tices. In this Report it is shown that modulation of the 
optical lattice along the longitudinal direction enables to 
realize the phenomenon of many-body CDT in a purely 
classical setting. As compared to quantum simulators 
of the Bose-Hubbard dynamics, in which the measurable 
quantities are usually the 'macroscopic' boson number 
and phase difference in the two wells [18], the proposed 
optical setting offers the rather unique possibility of a 
direct visualization of the full many-body tunneling dy- 



namics in Fock space. Moreover, it offers the ability to 
prepare and observe small systems of particle numbers, 
allowing to study the transition between few-body and 
many-body behavior. 

The starting point of our analysis is provided by a 
standard two-mode Bose-Hubbard Hamiltonian, which 
describes the tunneling dynamics of bosons in a double- 
well potential (see, for instance, [9t llSj) 

H{t) = Y(aja^ -f a\ai) -f ^^(a|a; - aXhrf , (1) 

where r (right) and / (left) arc the well sites, Ofe and a|, 
(k — l,r) are the bosonic annihilation and creation op- 
erators, V describes the constant tunneling rate between 
the two modes, and g(t) is the modulated interaction 
strength between same-mode bosons. The total number 
of bosons N = ajai + dja^ is a conserved quantity and 
the dimension of the Hilbert space is {N + 1). A sinu- 
soidal modulation g(t) = gisin^ut), without dc bias, is 
assumed for the interaction strength [19]. To study the 
phenomenon of many-body CDT, in Ref.[^ a Floquet 
analysis of the Hamiltonian H(t) was performed after in- 
troduction of the Schwinger representation of the angular 
momentum operators. To provide an optical realization 
of the two-mode Bose-Hubbard Hamiltonian and to un- 
derstand the many-particle CDT control scheme in an 
optical tunneling setup, we consider here the Fock space 
representation of the Hamiltonian H(t). If we expand 
the vector state of the system IV'(O) on the basis of Fock 
states with constant particle number N, i.e. after setting 

the evolution equations for the amplitude probabilities 
ci{t) to find / bosons in the left well and the other {N — l) 
bosons in the right well, as obtained from the Schrodinger 
equation ihdt\4'{t)) = H\tp{t)), read explicitly [2D] 

dci 

i-rr = + Ki^lCi-i) + g{t)ViCi (3) 

at 

{I — 0, 1, 2, N), where we have set 

ni = lVil + i)iN~l), Vi^li2l-Nr. (4) 
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The normalization condition X]/^o = 1 holds. 

The evolution equations (3) for the Fock state am- 
plitudes ci can be viewed as formally analogous to 
the coupled-mode equations describing light transport 
in a tight-binding array composed by (N + 1) waveg- 
uides with longitudinally-modulated propagation con- 
stant shift g{t)Vi and engineered coupling rates k; be- 
tween adjacent waveguides, in which the temporal evo- 
lution of the Fock-state amplitudes of the Bose-Hubbard 
Hamiltonian is mapped into the spatial evolution of the 
modal field amplitudes of light waves in the various 
waveguides along the axial direction t [T7]. The frac- 
tional light power distribution |c;p in the various waveg- 
uides of the array at the propagation distance t thus re- 
produces the distribution of occupation probabilities of 
the bosons between the two wells of the bistable poten- 
tial. It is worth noticing that, in the absence of par- 
ticle interaction, i.e. for g — 0, the lattice model (3) 
associated to the Bose-Hubbard Hamiltonian was previ- 
ously introduced in the photonic context to realize ex- 
act spatial beam self-imaging in finite waveguide arrays 
[3T] and shown to belong to a rather general class of 
exactly-solvable self-imaging tight-binding lattices with 
equally-spaced energy levels [25]. In the high-frequency 
modulation regime, the many-body CDT control method 
suggested in Ref. [5^ can be simply explained as a tunnel- 
ing inhibition process between certain waveguides in the 
array. After introduction of the slowly-varying ampli- 
tudes ai = ci exp[— z application of the averag- 
ing method yields the following evolution equations 
for the amplitudes ai 



. dai 

'' dt 



where we have set 

ai = Ki (exp 







(5) 



(6) 



and the bracket stands for the time average over one 
oscillation cycle. Equations (5) and (6) provide a first- 
order approximation to the exact tunneling dynamics in 
the lattice (for more details see [23]). For a sinusoidal 
modulation g{t) = giam{ujt) and using Eq.(4), the ex- 
plicit form of the effective coupling rate <ti reads 



ai = KieTq){iipi)Jo 



giiN~2l-l) 



(7) 



where Jq is the Bessel function of order zero and tpi = 
{gi/ijj){2l — N — 1) — Itt/2 is an inessential phase term. 
Suppose now that the system is initially prepared with 
all the N bosons on the right well, i.e. q — Si^q. Ac- 
cording to Ref.[n|, a desired number Iq — 0,1,2,... of 
particles can be allowed to tunnel to the left well pro- 
vided that the amplitude gi of the modulation is chosen 
in such a way that gi{N — 2lo — 1)/uj is a root of the Jq 
Bessel function. In this case, since the effective hopping 
rate aig between the sites Iq and Iq + I vanishes, one has 
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FIG. 1. (Color online) Evolution of the occupation probabil- 
ities |c,(t)p of iV = 10 bosons for the two-site Bose-Hubbard 
Hamiltonian (1) for v — 0.08 mm~^, lj — 0.628 mm^^ 
and for a few values of the modulation amplitude gi: (a) 



gi = 1.679 X 10""* mm"^; (b) gi = 2.159 x 10"" mm"\ (c) 
gi = 3.022 X lO""* mm-\ and (d) gi = 1.346 x 10""* mm~\ 
The four values of the modulation amplitude gi correspond 
to the four vertical lines I-IV in the quasi-energy diagram of 
Fig. 2(b). The modulation amplitudes in (a), (b) and (c) cor- 
respond to selective tunneling of zero, one and two bosons, 
respectively, as schematically indicated in the panels at the 
top of the figures. 
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FIG. 2. (Color online) (a) Evolution of the normalized popu- 
lation imbalance S{t) corresponding to the simulations shown 
in Fig.l (curve 1: gi = 1.679 x lO""* mm~^; curve 2; 
gi = 2.159 X 10"* mm"^; curve 3: gi = 3.022 x 10"'' mm"^; 
curve 4: gi = 1.346 x 10"" mm"^). (b) Quasi-energy spec- 
trum of the two-site Bose-Hubbard Hamiltonian (1) versus gi 
for V = 0.08 mm"^ and cj — 0.628 mm~^. The vertical lines 
I, II and III correspond to selective tunneling of zero, one and 
two bosons, respectively, as shown in Figs. 1(a), (b) and (c). 



ci{t) = for / > + 1, and tunneling to the lattice sites 
I = lo + l,lo + 2,...,N is thus forbidden. Since |ci(i)p is 
the probability to find I bosons in the left well and the 
other {N — I) bosons in the right well, this means that 
no more than Iq particles can tunnel from the right to 
the left well. As discussed in Ref.|9|, the conditions for 
selective CDT correspond to the crossing or touching of 
quasi-energy Floquet states of Eq. (3) with opposite par- 
ity. In our optical realization of the driven Bose-Hubbard 
Hamiltonian, the onset of many-body CDT is thus sim- 
ply visualized as the inhibition of light tunneling in a 
portion of the array. As an example, Figs.l(a-c) show 
the evolution of the occupation probabilities |ci(i)p, nu- 
merically computed from Eq.(3), corresponding to tun- 
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neling of zero, one and two particles for a iV = 10 
bosonic system, initially prepared with all the particles 
in the right well, for parameter values v — 0.08 mm~^, 
oj = 0.628 mm~^ and for the three different values gi of 
the modulation satisfying the selective CDT conditions 
gi{N - 2lo -l)/uj = 2.405 with /q = 0, 1 and 2. The cor- 
responding evolution of the normalized population im- 

balance ^J, given by 8(1) = E^IoK^ " 20/^^] |c/(i)P, is 
depicted in Fig. 2(a). The selective CDT realized in the 
three cases corresponds to the three points I, II and III in 
the quasi-energy diagram shown in Fig. 2(b). For compar- 
ison, in Fig. 1(d) and curve 4 of Fig. 2(a) the behavior of 
the occupation probabilities and normalized population 
imbalance is also shown for a wrong value of the mod- 
ulation amplitude gi, which does not correspond to de- 
generacy of Floquet states [see hne IV in Fig. 2(b)]. Note 
that, as opposed to Figs.l(a)-(c), in this case the light 
waves spread over the entire waveguides of the array. 

To realize the lattice model (3), the hopping rates and 
site potentials should be tailored according to Eqs.(4). In 
a waveguide array, this requires a control of waveguide 
spacing and of refractive index contrasts and/or channel 
sizes of the various waveguides. The tight-binding model 
(3) can be derived using a variational procedure starting 
from the paraxial and scalar wave equation for the elec- 
tric field amplitude 4){x,t) describing the propagation of 
monochromatic light waves at wavelength A in an array 
of (A'' -t- 1) waveguides with a longitudinally- modulated 
refractive index profile n{x, t) and substrate index Us 



iXdtct) = -^dl(p+[ns - n{x,t)](l), 
2n., 



(8) 



where A = A/(27r) is the reduced wavelength of pho- 
tons and t is the longitudinal spatial coordinate. The 
hopping rates fc; and propagation constant shifts gVi 
entering in the reduced equations (3) are given by 
certain overlapping integrals involving the modal field 
and index profiles of the guides (for details see, for 
instance, j24j ) . In particular, to independently engineer 
the coupling constants ki and propagation constant 
shifts gVi of the waveguides, one can assume a chain 
of waveguides with equal normalized refractive index 
profile G{x), but with different and longitudinally- 
modulated index contrasts Ani{t) {I = 0,1,2, ...,N) 
and spacing di — xi — xi^i {I = l,2,...,N), i.e. 
n{x, t) — Us = X^i^o A'n'i{t)G{x — Xi). To realize the two- 
mode Bose-Hubbard lattice (3) , the values of di and An; 
slightly vary around some mean values dr and An that 
define a uniform array ITi . The waveguide separation di 
mainly determines the value of the coupling rate 
with a characteristic exponential dependence of 
from di, whereas the index change An; mainly defines 
the propagation constant mismatch giVi. Numerical 
simulations were performed to reproduce the parame- 
ter conditions of Fig.l for an operational wavelength 
A = 633 and assuming n^ = 1.45, which applies to fused 
silica waveguide arrays. The normalized waveguide 
profile used in the simulations is given by G{x) = 




FIG. 3. (Color online) (a-d) Evolution of beam intensity 
(snapshot of \4>{x,t)\^) in four 10-cm-long longitudinally- 
modulated waveguide arrays, composed by (A'^ -|- 1) = 11 
waveguides, that realize the lattice models of Figs. 1 (a-d). 
The refractive index profiles n{x) — Us of the arrays, corre- 
sponding to the maximum/ minimum of modulation ampli- 
tude [sm(uiz) = ±1] are shown in the upper panels of the 
figures by the solid/dashed curves. The behavior of the nor- 
malized population imbalance S{t), retrieved from the beam 
center of mass position, is shown in (e) for the four arrays. 



{erf [(a: + w)/D,] - erf[(x - w)/D,]}/[2eTi{w/D,)], 
with channel width 2w — 4 /im and diffusion length 
Dx — 0.3 ^m~^. To determine the distribution of 
distances di, a reference value An = 2 x 10^"^ of 
refractive index contrast was assumed, and corre- 
spondingly the coupling rate k between two adjacent 
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waveguides versus distance d was computed, yielding 
to a good approximation the exponential dependence 
K{d) — Koexp[— 7(0? — dr)] for distances not too 
far from the reference distance = 9 /^m, where 
kq — 0.2144 mm~^ and 7 ~ 0.6 /im~^. The resulting 
distance distribution that realizes the coupling rates 
(4) with V = 0.08 mm~^ indicates that di varies in the 
range 9-9.916 /ini. By modulating the index contrasts 
An/ of the waveguides around the reference value An, 
four different arrays were then designed to realize the 
four different interaction regimes of Fig.l. A numerical 
computation of the propagation constant of the channel 
waveguide indicates that a propagation constant shift 
gi sm{ujt)Vi is approximately obtained by assuming a 
refractive index contrast An; = An + /Sgi sin(a;t) V/A, 
where An = 2 x 10"'^ is the reference value of the 
index contrast and /? is a numerical factor of order 
one {P ~ 1.23 for the chosen error function waveguide 
profile). The resulting refractive index profile of the 
arrays, corresponding to the maximum [sin(a;t) — 1] and 
minimum [sin(a;t) — —1] of the modulation cycle, are 
depicted in the upper panels of Figs.3(a-d). Such arrays 
could be fabricated in fused silica by the femtosecond 
laser writing technique |25| , in which the different refrac- 
tive index contrasts are obtained by varying the speed 
of the writing laser beam. Figures 3(a-d) show the evo- 
lution of light intensity |(/)(a;,i)p along the four arrays, 



as obtained by numerical simulations of Eq.(8) using a 
standard pseudospectral split-step method, when the left 
boundary waveguide is excited at the input plane. The 
behavior of the normalized population imbalance S{t)^ 
which provides a clear signature of selective CDT [9], 
can be retrieved from a measurement of the beam center 
of mass position {x{t)) = J dxx\(j){x,t)\'^ / J dx\(j){x,t)\'^ 
using the simple relation S{t) ~ 1 - 2{x{t)) / {N dr) [17] . 
The corresponding behavior of S{t) for the four arrays is 
depicted in Fig. 3(e), which clearly indicates the onset of 
selective CDT according to the results shown in Figs.l 
and 2. 

In conclusion, a classic wave optics realization of the 
selective many-body coherent destruction of tunneling, 
recently predicted for interacting many-boson systems 
in Ref.|9j, has been proposed for light transport in mod- 
ulated waveguide lattices. In such photonic structures 
light propagation in the various waveguides maps the 
evolution of the bosonic distribution in the double-well 
potential, and thus enables to visualize the entire 
many-body tunneling dynamics in the Fock space, a 
possibility which would be hard to achieve in matter 
wave systems. 
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